home *** CD-ROM | disk | FTP | other *** search
/ C/C++ Users Group Library 1997 August / Walnut Creek CDROM.7z / VOL_400 / 446_01 / DOC / DPBASICS / ODE / ODE.C < prev    next >
Encoding:
C/C++ Source or Header  |  1996-04-18  |  1.9 KB  |  67 lines
  1. #include <ode.h>
  2. ODESolver::  ODESolver (ODESolverUDC& f)  
  3. : eqdef(f) // f, or eqdef as in <ode.h>, is just a dummy name
  4. { init(); }  // initializations performed in init
  5.  
  6. void ODESolver:: init ()
  7. { scratch1.redim (eqdef.size()); }
  8.  
  9. ForwardEulerODE:: ForwardEulerODE (ODESolverUDC& eqdef)
  10. : ODESolver (eqdef) {}
  11.  
  12. void ForwardEulerODE:: advance (Vec(real)& y, real& t, real& dt)
  13. {
  14.   eqdef.equation (scratch1, y, t);
  15.   for (int i = 1; i <= y.size(); i++)
  16.     y(i) += dt * scratch1(i);
  17.   t += dt;
  18. }
  19.  
  20. RungeKutta4ODE:: RungeKutta4ODE (ODESolverUDC& eqdef)
  21. : ODESolver (eqdef), scratch2(eqdef.size()), scratch3(eqdef.size()) {}
  22.  
  23. void RungeKutta4ODE:: advance (Vec(real)& y, real& t, real& dt) 
  24. {
  25.   real dt2 = 0.5*dt;
  26.   real dt6 = dt/6.0;
  27.   eqdef.equation (scratch1, y, t);
  28.   for (int i = 1; i <= y.size(); i++)
  29.     scratch2(i) = y(i)  +  dt2 * scratch1(i);
  30.   eqdef.equation (scratch1, scratch2, t+dt2);
  31.   for (i = 1; i <= y.size(); i++)
  32.     scratch2(i) = y(i)  +  dt2 * scratch1(i);
  33.   eqdef.equation (scratch3, scratch2, t+dt2);
  34.   for (i = 1; i <= y.size(); i++)
  35.     {
  36.       scratch2(i) = y(i)   +  dt * scratch3(i);
  37.       scratch3(i) = scratch1(i) + scratch3(i);
  38.     }
  39.   eqdef.equation (scratch1, scratch2, t+dt);
  40.   eqdef.equation (scratch2, y, t); 
  41.  
  42.   for (i = 1; i <= y.size(); i++)
  43.     y(i) = y(i)  +  dt6 * ( scratch1(i) + scratch2(i) + 2*scratch3(i) );
  44.   t += dt;
  45. }
  46.  
  47. void Oscillator:: equation (Vec(real)& f, const Vec(real)& y, real t)
  48. {
  49.   f(1) = y(2);
  50.   f(2) = -c1*(y(2)+c2*y(2)*abs(y(2)))-c3*(y(1)+c4*pow3(y(1))) + sin (omega*t);
  51. }
  52.  
  53. void Oscillator:: scan (Is is)
  54. {
  55.   // reades c1, c2, c3, c4 and omega:
  56.   is >> c1 >> c2 >> c3 >> c4 >> omega;
  57.   s_o << "\nHere are the Oscillator parameters:\n c1=" << c1 << " c2=";
  58.   s_o << c2 << " c3=" << c3 << " c4=" << c4 << " omega=" << omega << '\n';
  59. }
  60.  
  61. void Oscillator:: print (Os os)
  62. {
  63.   os << aform("d^2y/dt^2 + %g*(dy/dt + %g*dy/dt*|dy/dt|) + %g*(y + %g*y^3) \
  64.         = sin(%g*t)\n",c1,c2,c3,c4,omega);
  65. }
  66.  
  67.